In [1]:
#PSD on profilometer of slumped parts

%reset
%load_ext autoreload

%autoreload 2
%matplotlib qt

Once deleted, variables cannot be recovered. Proceed (y/[n])? y
In [ ]:
#plot psd of calibration flat to get signature of instrument
#compare profilometer data of OAB103 different measurement with different coatings
#add CCI OAB mandrel data from first set
#WFS data on profilometer and mid-frequencies
#compare psd of slumping models at different temperatures

%qtconsole
In [2]:
import matplotlib.pyplot as plt
import numpy as np
import os
from astropy.io import fits
from dataIO.span import span
from pyGeneralRoutines.fn_add_subfix import fn_add_subfix

from IPython.display import display

#from skimage import data as skidata
#from skimage import img_as_float

#from plot_grains import *
#from radialProfiles import *
#from peak_finder import *

from pyProfile.psd import plot_sig_psd
from psd2d import *

from pySurf.points import *

## trick to import from a non released library
import sys
sys.path.append(r'..\OP1surface')
from fit_cylinder import *

outfolder=r'PSDtest'

#this is the simplest way to make ipython print numbers in non scientific format 
%precision '%f'  
np.set_printoptions(formatter={'float_kind':lambda x: '%f'%x})
In [103]:
def correctedPSD(pts,crop=None,step=None,nsigma=None,outname=None):
    """Perform a set of operations on pts and plot intermediate 
    steps:
    rebin based on step, this must be provided to convert from pointcloud to matrix
    remove outliers line by line on the base of sigma
    if outname is provided, save figures.
    Return 1d frequency and matrices of psds along axis before and after correction.
    """
    if outname is not None:
        outfolder=os.path.dirname(outname)
        outname=os.path.basename(outname)
    
    plt.figure(1)
    plt.clf()
    plot_points(pts,scatter=1,aspect='auto')
    display(plt.gcf())

    if crop is not None:
        pts=crop_points(pts,xrange=crop[0],yrange=crop[1])
        
    pts=level_points(pts)
    

    print "rebin with step %s,%s"%step
    # rebin
    ## create a grid for bins
    s=span(pts[:,0:2],axis=0)
    xb=np.arange(s[0][0],s[0][1],step[0])  #x and y bins
    yb=np.arange(s[1][0],s[1][1],step[1])
    p,x,y=rebin_points(pts,bins=(xb,yb),matrix=True)  

    plt.figure(2)
    plt.clf()
    plt.imshow(p,extent=np.hstack([span(x),span(y)]),aspect='equal',
               vmin=-10,interpolation='none',
                origin='lower')
    plt.colorbar()
    if outname is not None:
        plt.savefig(os.path.join(outfolder,fn_add_subfix(outname,
                                                         '_data','.png')))
    display(plt.gcf())

    if nsigma is not None:
        #remove outliers and plot comparison
        pl=remove_outliers2d(p,fignum=3,nsigma=nsigma)
        if outname is not None:
            plt.savefig(os.path.join(outfolder,fn_add_subfix(outname,'_datacorr','.png')))
        display(plt.gcf())
    else:
        pl=p
    
    f,pp1,pp2=compare_2dpsd(p,pl,x,y,fignum=5)
    if outname is not None:
        plt.savefig(os.path.join(outfolder,fn_add_subfix(outname,'_2dPSDcorr','.png')))
    display(plt.gcf())
    
    ap1=avgpsd2d(pp1)
    ap2=avgpsd2d(pp2)

    plt.figure(4)
    plt.clf()
    #plt.title('Calibration mandrel')
    plt.plot(f[1:],ap1[1:-1],'r',label='Raw')
    plt.plot(f[1:],ap2[1:-1],'b',label='Outlier removed')
    plt.loglog()
    plt.ylabel('um^2')
    plt.xlabel('mm^-1')
    plt.legend(loc=0)
    plt.grid(1)
    if outname is not None:
        plt.savefig(os.path.join(outfolder,fn_add_subfix(outname,'_PSD','.png')))
    display(plt.gcf())
    return f,ap1,ap2

Calibration mandrel

In [104]:
#try to repeat results with a function
df=r'../OP1surface/measure_data/calibration_m3/2013_12_23/OAB150x150_01_Height.txt'
pts=get_points(df,center=(0,0),scale=(-1,1.,1))

crop=[[-60,60],[-60,60]]
step=(0.5,0.5)
nsigma=1.5
#pts=get_points(df,center=(0,0),scale=(1,1,1),delimiter=' ')
f,ap1,ap2=correctedPSD(pts,crop=crop,step=step,nsigma=nsigma)

plt.savefig(os.path.join(outfolder,fn_add_subfix(os.path.basename(df),'_PSD','.png')))
np.savetxt(os.path.join(outfolder,fn_add_subfix(os.path.basename(df),'_PSD','.dat')),
           np.vstack([f[:],ap1[:-1],ap2[:-1]]).T)

rebin with step 0.5,0.5

OAB mandrel Processed in same way.

In [99]:
#try to repeat results with a function
df=r'C:\Users\Vincenzo\Google Drive\POOL\OP1surface\OAB1\xysurface150x150_Height_transformed_deltaR_4inch.dat'


crop=None
step=(0.5,0.5)
nsigma=1.5

pts=get_points(df,center=(0,0),scale=(1,1,1),delimiter=' ')
f,ap1,ap2=correctedPSD(pts,crop=crop,step=step,nsigma=nsigma)

plt.savefig(os.path.join(outfolder,fn_add_subfix(os.path.basename(df),'_PSD','.png')))
np.savetxt(os.path.join(outfolder,fn_add_subfix(os.path.basename(df),'_PSD','.dat')),
           np.vstack([f[:],ap1[:-1],ap2[:-1]]).T)

rebin with step 0.5,0.5

Mandrel OP1

In [33]:
#Here I start from raw data, need to fit cylinder

df=r'..\OP1surface\OP1\09_op1_yxsurf_Height_transformed.dat'
#measure_data\2015_07_10\06_OP1_xyscan_Height.txt'

crop=[[-23,23],[-50,50]]
step=(0.5,0.5)
nsigma=1.5

pts=get_points(df,center=(0,0),scale=(1,1,1),delimiter=' ') #z-1 per mandreino
pts=crop_points(pts,crop[0],crop[1])
In [34]:
plt.clf()
a=plot_points(pts,scatter=True,aspect='equal')
plt.grid(1)
display(plt.gcf())
In [35]:
#here i can remove cylinder
# wrapper around minimize.

def fit(pts,odr2,fit_func):
    """   fit pts and return residuals. Info are printed. odr2 is starting guess.
    typically:
    fit_func=cylinder_error3 #function to be minimized
    odr2=(span(xg).sum()/2.,220.,0,0.) #use nominal value for guess direction
    deltaR=fit_func(pts,odr2,fit_func).
    pts and deltaR enter and exit with z in microns
    """
    
    pts=pts/[1,1,1000.] #data are um, convert to mm before fit
    result=minimize(fit_func,x0=(odr2),
                    args=(pts,),
                    options={'maxiter':500},
                    method='Nelder-Mead')#,callback=p)
    print '-----------------------------------'
    print 'Results of fit on region:'
    print result    
    odr=result.x
    fom,deltaR,coeff=fit_func(odr,pts,extra=True)
    deltaR[:,2]=deltaR[:,2]*1000
    print 'Angle of cyl axis with y axis (deg):'
    print np.arccos(np.sqrt(1-(odr[-2:]**2).sum()))*180/np.pi
    print 'Radius:'
    print coeff
    print 'rms residuals=%s um'%(fom*1e3)
    return deltaR
In [39]:
odr2=(span(pts[:,0]).sum()/2.,220.,0,0.) #use nominal value for guess direction
deltaR=fit(pts,odr2,fit_func=cylinder_error3) #cylinder_error3 #function to be minimized

-----------------------------------
Results of fit on region:
 final_simplex: (array([[0.000000, 220.457499, 0.004209, -0.000006],
       [0.000000, 220.457548, 0.004210, -0.000006],
       [0.000000, 220.457590, 0.004209, -0.000006],
       [0.000000, 220.457588, 0.004209, -0.000006],
       [0.000000, 220.457430, 0.004209, -0.000006]]), array([0.003944, 0.003944, 0.003944, 0.003944, 0.003944]))
           fun: 0.0039443651906279483
       message: 'Optimization terminated successfully.'
          nfev: 168
           nit: 96
        status: 0
       success: True
             x: array([0.000000, 220.457499, 0.004209, -0.000006])
Angle of cyl axis with y axis (deg):
0.241184635264
Radius:
221.378669474
rms residuals=3.94436519063 um
In [40]:
sigma=np.std(deltaR[:,2])
print sigma
plt.clf()
plot_points(deltaR,scatter=True,aspect='equal',vmin=-sigma,vmax=sigma)
display(plt.gcf())
plt.grid(1)

3.94436519063
In [41]:
#remove outliers based on deltaR from pts

sigma=np.std(deltaR[:,2])
mean=np.mean(deltaR[:,2])
mask=(deltaR[:,2]-mean)**2<(nsigma*sigma)**2 #points to keep
pts=pts[mask,:]
print sigma,np.std(deltaR[mask,2])

3.94436519063 0.754320797265
In [43]:
odr2=(span(pts[:,0]).sum()/2.,220.,0,0.) #use nominal value for guess direction
deltaR=fit(pts,odr2,fit_func=cylinder_error3) #cylinder_error3 #function to be minimized

-----------------------------------
Results of fit on region:
 final_simplex: (array([[0.000000, 220.455726, 0.004207, -0.000006],
       [0.000000, 220.455776, 0.004207, -0.000006],
       [0.000000, 220.455663, 0.004207, -0.000006],
       [0.000000, 220.455642, 0.004207, -0.000006],
       [0.000000, 220.455637, 0.004207, -0.000006]]), array([0.000754, 0.000754, 0.000754, 0.000754, 0.000754]))
           fun: 0.00075429435148178169
       message: 'Optimization terminated successfully.'
          nfev: 169
           nit: 92
        status: 0
       success: True
             x: array([0.000000, 220.455726, 0.004207, -0.000006])
Angle of cyl axis with y axis (deg):
0.241055005866
Radius:
221.376888299
rms residuals=0.754294351482 um
In [44]:
sigma=np.std(deltaR[:,2])
print sigma
plt.clf()
plot_points(deltaR,scatter=True,aspect='equal',vmin=-sigma,vmax=sigma)
display(plt.gcf())
plt.grid(1)

0.754294351482
In [45]:
f,ap1,ap2=correctedPSD(deltaR,crop=crop,step=step,nsigma=nsigma,outname=df)
In [46]:
plt.savefig(os.path.join(outfolder,fn_add_subfix(os.path.basename(df),'_PSD','.png')))
np.savetxt(os.path.join(outfolder,fn_add_subfix(os.path.basename(df),'_PSD','.dat')),
           np.vstack([f[:],ap1[:-1],ap2[:-1]]).T)

Plot all together

In [108]:
df1=r'../OP1surface/measure_data/calibration_m3/2013_12_23/OAB150x150_01_Height.txt'
df2=r'C:\Users\Vincenzo\Google Drive\POOL\OP1surface\OAB1\xysurface150x150_Height_transformed_deltaR_4inch.dat'
df3=r'..\OP1surface\OP1\09_op1_yxsurf_Height_transformed.dat'

f1,p1,p1c=np.genfromtxt(os.path.join(outfolder,fn_add_subfix(
            os.path.basename(df1),'_PSD','.dat')),
            unpack=True)
f2,p2,p2c=np.genfromtxt(os.path.join(outfolder,fn_add_subfix(
            os.path.basename(df2),'_PSD','.dat')),
            unpack=True)
f3,p3,p3c=np.genfromtxt(os.path.join(outfolder,fn_add_subfix(
            os.path.basename(df3),'_PSD','.dat')),
            unpack=True)

plt.figure('PSD')
plt.clf()
plt.title('PSD of mandrels')
plt.plot(f1[1:],p1[1:],'b--',label='Ref. flat Raw')
plt.plot(f1[1:],p1c[1:],'b',label='Ref. flat outlier removed')
plt.plot(f2[1:],p2[1:],'r--',label='OAB1 Raw')
plt.plot(f2[1:],p2c[1:],'r',label='OAB1 outlier removed')
plt.plot(f3[1:],p3[1:],'g--',label='OP1 Raw')
plt.plot(f3[1:],p3c[1:],'g',label='OP1 outlier removed')
plt.loglog()
plt.ylabel('um^2')
plt.xlabel('mm^-1')
plt.xlim((0.8e-2,2e0))
plt.legend(loc=0)
plt.grid(1)
display(plt.gcf())

plt.savefig(os.path.join(outfolder,'mandrels_PSD.png'))
In [109]:
#corrected only
plt.figure('PSD')
plt.clf()
plt.title('PSD of mandrels')
#plt.plot(f1[1:],p1[1:],'b--',label='Ref. flat Raw')
plt.plot(f1[1:],p1c[1:],'b',label='Ref. flat')
#plt.plot(f2[1:],p2[1:],'r--',label='OAB1 Raw')
plt.plot(f2[1:],p2c[1:],'r',label='OAB1')
#plt.plot(f3[1:],p3[1:],'g--',label='OP1 Raw')
plt.plot(f3[1:],p3c[1:],'g',label='OP1')
plt.loglog()
plt.ylabel('um^2')
plt.xlabel('mm^-1')
plt.xlim((0.8e-2,2e0))
plt.legend(loc=0)
plt.grid(1)
display(plt.gcf())

plt.savefig(os.path.join(outfolder,'mandrels_PSD_corronly.png'))